#Author : Floriane Simonet
#e-mail : simonet.floriane@gmail.com
#! usr/bin/perl -w
use strict;



#sub create_snp_set2{

#	my($ref_listchrp,$ref_snpfile,$nb_gp,$setfile1)=@_;
#	my @listchrp=@$ref_listchrp;
#	my %snpfile=%$ref_snpfile;
#	my $currentgp=0;
#	my $sumPos=0;
#	my $cpt=0;
#	my $meanPos=0;
#	my $posm1=0;
#	my %tabpos=();
	

	# write newbimfile to impute
#	my $setfile2="snp_setfile2.txt";
#	open(SETFILE2,">$setfile2") or die "cannot open SETFILE2 : $setfile2 $!\n";
#	open(SETFILE1,$setfile1) or die "cannot open SETFILE1 : $setfile1 $!\n";
	
#	while(my $line=<SETFILE1>){

#		chomp($line);
#		my($chr,$pos,$rs,$gp)=split(/\t/,$line);
		
#		print "$chr\t $pos \t $rs \t $gp \n";
#		if(($gp!=$currentgp) && ($cpt!=0)){ #new group
#			print "2\n";
#			# save information about last group
#			print "currrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrent : $currentgp\n";
#			$tabpos{$currentgp}{'lastPos'}=$posm1;
#			$meanPos= sprintf("%.0f", ($sumPos/$cpt));
#			print "******************$meanPos\t$sumPos\t$cpt\t$posm1\n";
#			$tabpos{$currentgp}{'meanPos'}=$meanPos;

			# start new group
#			$tabpos{$gp}{'firstPos'}=$pos;
#			$sumPos=$pos;
#			$cpt=1;
#			$currentgp=$gp;
#		}
#		elsif(($gp==$currentgp) && ($gp!=0)){ #same group
#			print "3\n";
#			$sumPos=$sumPos+$pos;
#			$cpt=$cpt+1;
#		}
#		else{ #first group
#			print "1\n";
#			$sumPos=$pos;
#			$tabpos{$gp}{'firstPos'}=$pos;
#			$cpt=1;
#			$currentgp=$gp;
#		}
#	$posm1=$pos;
#	print "\n meanpos=> $tabpos{$currentgp}{'meanPos'}\n";
#	print "\n firstpos=> $tabpos{$gp}{'firstPos'}\n";
#	print "\n lastpos=> $tabpos{$currentgp}{'lastPos'}\n";
#	}
#	close SETFILE1;
#	close SETFILE2;
#	print "$tabpos{2}{'meanPos'}\n";
#	print "$tabpos{2}{'firstPos'}\n";
#	print "$tabpos{2}{'lastPos'}\n";
#}


#	for(my $i=0; $i < $nb_gp; $i++){
#		}
#	for(my $i=0; $i < $nb_gp; $i++){ 
#		my $chp=$listchrp[$i]; 
#	}

# my $mean = 
# my $roundmean = sprintf("%.0f", $mean); => pour faire un arrondi
#		print $snpfile{$chp}{'chr'}."\n";
#		print $snpfile{$chp}{'rsID'}."\n";
#		print $snpfile{$chp}{'cMpos'}."\n";
#		print $snpfile{$chp}{'BPpos'}."\n";
#		print $snpfile{$chp}{'a1'}."\n";
#		print $snpfile{$chp}{'a2'}."\n";
#		if((exists($snpfile{$listchrp[$j-1]})) && ($snpfile{$chp}{'chr'} == $snpfile{$listchrp[$j-1]}{'chr'})){ # same chr
#			$D=$pos1 - $snpfile{$chp}{'BPpos'};
#			if( ((($snpfile{$listchrp[$j]}{'BPpos'}) - ($snpfile{$listchrp[$j-1]}{'BPpos'})) < 1000000) && ($D < 4000000) ){ # si la distance du marqueur avec celui d'avant est < 1MB et si la distane du premier marqueur du groupe à ce marqueur ci est < 4Mb
#				print SETFILE2 "$snpfile{$chp}{'chr'}\t$snpfile{$chp}{'BPpos'}\t$snpfile{$chp}{'rsID'}\t$gp\n";
#			}
#			else{
#				$gp++;
#				print SETFILE2 "$snpfile{$chp}{'chr'}\t$snpfile{$chp}{'BPpos'}\t$snpfile{$chp}{'rsID'}\t$gp\n";
#				$pos1=$snpfile{$chp}{'BPpos'};
#			}
#		}
#		elsif ((exists($snpfile{$listchrp[$j-1]})) && ($snpfile{$chp}{'chr'} != $snpfile{$listchrp[$j-1]}{'chr'} ) || (($j-1) == -1)){ # new chr
#			$gp++;
#			print SETFILE2 "$snpfile{$chp}{'chr'}\t$snpfile{$chp}{'BPpos'}\t$snpfile{$chp}{'rsID'}\t$gp\n";
#			$D=$snpfile{$chp}{'BPpos'};
#		}	


sub create_snp_set2{

	my($ref_listchrp,$ref_snpfile,$nb_gp,$ref_snpset1)=@_;
	my @listchrp=@$ref_listchrp;
	my %snpfile=%$ref_snpfile;
	my %setfile1=%$ref_snpset1;


	foreach my $rsid (keys(%setfile1)){
				print "#########################\n";
				print "chr : $setfile1{$rsid}{'chr'}\n";
				print "pos1 : $setfile1{$rsid}{'pos1'}\n";
				print "pos2 : $setfile1{$rsid}{'pos2'}\n";
				print "gp : $setfile1{$rsid}{'gp'}\n";
	}
	
	close SETFILE1;
	close SETFILE2;
}

##########

sub create_snp_set1{

	my($ref_listchrp,$ref_snpfile)=@_;
	my @listchrp=@$ref_listchrp;
	my %snpfile=%$ref_snpfile;

	my%snpset1=();
	my $gp=0;
	my $D=0;
	my $pos1=0;
	my $tmprsid;
#	my $d=0;

	# write newbimfile to impute
	my $setfile1="snp_selected_setfile.txt";
	open(SETFILE1,">$setfile1") or die "cannot open SETFILE1 : $setfile1 $!\n";

	for(my $j=0; $j < scalar(@listchrp); $j++){ 

		my $chp=$listchrp[$j]; 
#		$d=(($snpfile{$listchrp[$j]}{'BPpos'}) - ($snpfile{$listchrp[$j-1]}{'BPpos'}));
#		print $snpfile{$chp}{'chr'}."\n";
#		print $snpfile{$chp}{'rsID'}."\n";
#		print $snpfile{$chp}{'cMpos'}."\n";
#		print $snpfile{$chp}{'BPpos'}."\n";
#		print $snpfile{$chp}{'a1'}."\n";
#		print $snpfile{$chp}{'a2'}."\n";
		if((exists($snpfile{$listchrp[$j-1]})) && ($snpfile{$chp}{'chr'} == $snpfile{$listchrp[$j-1]}{'chr'})){ # same chr
			$D=$pos1 - $snpfile{$chp}{'BPpos'};
			if( ((($snpfile{$listchrp[$j]}{'BPpos'}) - ($snpfile{$listchrp[$j-1]}{'BPpos'})) < 1000000) && ($D < 4000000) ){ # si la distance du marqueur avec celui d'avant est < 1MB et si la distance du premier marqueur du groupe à ce marqueur est < 4Mb
				print SETFILE1 "$snpfile{$chp}{'chr'}\t$snpfile{$chp}{'BPpos'}\t$snpfile{$chp}{'rsID'}\t$gp\n";
#				print "###############\n$snpfile{$chp}{'chr'}\t$snpfile{$chp}{'BPpos'}\t$snpfile{$chp}{'rsID'}\t$gp\n";
			}
			else{
				$snpset1{$tmprsid}{'pos2'}=$snpfile{$listchrp[$j-1]}{'BPpos'} if(exists($snpfile{$listchrp[$j-1]}));
#				print " pos2 : $snpset1{$tmprsid}{'pos2'}	$tmprsid\n";
				$gp++;
				print SETFILE1 "$snpfile{$chp}{'chr'}\t$snpfile{$chp}{'BPpos'}\t$snpfile{$chp}{'rsID'}\t$gp\n";
#				print "###############\n$snpfile{$chp}{'chr'}\t$snpfile{$chp}{'BPpos'}\t$snpfile{$chp}{'rsID'}\t$gp\n";
				$snpset1{$snpfile{$chp}{'rsID'}}{'pos1'}=$snpfile{$chp}{'BPpos'};
				$snpset1{$snpfile{$chp}{'rsID'}}{'chr'}=$snpfile{$chp}{'chr'};
				$snpset1{$snpfile{$chp}{'rsID'}}{'gp'}=$gp;
				$tmprsid=$snpfile{$chp}{'rsID'};
#				print " pos1 : $snpset1{$snpfile{$chp}{'rsID'}}{'pos1'}	$snpfile{$chp}{'rsID'}\n";
#				print " chr : $snpset1{$snpfile{$chp}{'rsID'}}{'chr'}	$snpfile{$chp}{'rsID'}\n";
#				print " gp : $snpset1{$snpfile{$chp}{'rsID'}}{'gp'}	$snpfile{$chp}{'rsID'}\n";
				$pos1=$snpfile{$chp}{'BPpos'};
			}
		}
		elsif ((exists($snpfile{$listchrp[$j-1]})) && ($snpfile{$chp}{'chr'} != $snpfile{$listchrp[$j-1]}{'chr'} ) || (($j-1) == -1)){ # new chr
#			$snpset1{$tmprsid}{'pos2'}=$snpfile{$listchrp[$j-1]}{'BPpos'} if(exists($snpfile{$listchrp[$j-1]}) && !defined($snpset1{$tmprsid}{'pos2'}));
			$gp++;
			print SETFILE1 "$snpfile{$chp}{'chr'}\t$snpfile{$chp}{'BPpos'}\t$snpfile{$chp}{'rsID'}\t$gp\n";
#			print "###############\n$snpfile{$chp}{'chr'}\t$snpfile{$chp}{'BPpos'}\t$snpfile{$chp}{'rsID'}\t$gp\n";
			$snpset1{$snpfile{$chp}{'rsID'}}{'pos1'}=$snpfile{$chp}{'BPpos'};
			$snpset1{$snpfile{$chp}{'rsID'}}{'chr'}=$snpfile{$chp}{'chr'};
			$snpset1{$snpfile{$chp}{'rsID'}}{'gp'}=$gp;
			$tmprsid=$snpfile{$chp}{'rsID'};
#			print " pos1 : $snpset1{$snpfile{$chp}{'rsID'}}{'pos1'}	$snpfile{$chp}{'rsID'}\n";
#			print " chr : $snpset1{$snpfile{$chp}{'rsID'}}{'chr'}	$snpfile{$chp}{'rsID'}\n";
#			print " gp : $snpset1{$snpfile{$chp}{'rsID'}}{'gp'}	$snpfile{$chp}{'rsID'}\n";
			$D=$snpfile{$chp}{'BPpos'};
		}	
	}
	return ($gp,\%snpset1);	
	close(SETFILE1);
}

##########

sub sort_by_chrpos{
	
	my $ref_chrpos=$_[0];
	my @chrpos=@$ref_chrpos;

	my @sortchrpos = sort {
	my @a = $a =~/(\d+)\-(\d+)/;
	my @b = $b =~/(\d+)\-(\d+)/; 
	$a[0] <=> $b[0] || $a[1] <=> $b[1] ;
	} @chrpos; # sort by chromosome then sort by position

	return(\@sortchrpos);
}

##########

sub read_SNPselectfile{

	my $SNPselectfile=$_[0];
	my %snpfile=();
	my @chrpos=();

	open(SNPSELECTFILE, $SNPselectfile) or die "cannot open $SNPselectfile";
	while(my $line=<SNPSELECTFILE>)
	{
		chomp($line);
		my($chr,$rsID,$cMpos,$BPpos,$a1,$a2)=split(/\t/,$line); # read snp file
		
		my $chp=$chr."-".$BPpos; # get chr-position
		push(@chrpos,$chp);

		$snpfile{$chp}{'chr'}=$chr;
		$snpfile{$chp}{'rsID'}=$rsID;
		$snpfile{$chp}{'cMpos'}=$cMpos;
		$snpfile{$chp}{'BPpos'}=$BPpos;
		$snpfile{$chp}{'a1'}=$a1;
		$snpfile{$chp}{'a2'}=$a2;

	}
	close SNPSELECTFILE;

	my $ref_listchrp=&sort_by_chrpos(\@chrpos);
	my @listchrp=@$ref_listchrp;

	my ($nb_gp,$ref_snpset1)=&create_snp_set1(\@listchrp,\%snpfile);

	return ($nb_gp,$ref_snpset1);
}

##########

sub read_SNPgenotypedfile{
	
	my $filename=$_[0];
	my $nb_gp=$_[1];
	my $ref_snpset1=$_[2];
	my %file=();
	my @chrpos=();

	open(FILE, $filename) or die "cannot open $filename";
	while(my $line=<FILE>)
	{
		chomp($line);
		my($chr,$rsID,$cMpos,$BPpos,$a1,$a2)=split(/\t/,$line); # read snp file
		
		my $chp=$chr."-".$BPpos; # get chr-position
		push(@chrpos,$chp);

		$file{$chp}{'chr'}=$chr;
		$file{$chp}{'rsID'}=$rsID;
		$file{$chp}{'cMpos'}=$cMpos;
		$file{$chp}{'BPpos'}=$BPpos;
		$file{$chp}{'a1'}=$a1;
		$file{$chp}{'a2'}=$a2;
		
	}
	close FILE;

	my $ref_listchrp=&sort_by_chrpos(\@chrpos);
	my @listchrp=@$ref_listchrp;

	&create_snp_set2(\@listchrp,\%file,$nb_gp,$ref_snpset1);
}

############################################################

sub main(){
	
	my $SNPselectfile=$ARGV[0]; # SNPs selected for replication (target SNP)
	my $SNPgenotypedfile=$ARGV[1];	

	if(-f $SNPselectfile){
		my ($nb_gp,$ref_snpset1)=&read_SNPselectfile($SNPselectfile);
		&read_SNPgenotypedfile($SNPgenotypedfile,$nb_gp,$ref_snpset1);
	}
	else{
			print "$SNPselectfile doesn't exist.\n";
	}
exit;
}

main();

# command line example : perl imputation.pl ../newbimfile_mitral.input_10252011.out
